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Abstract 

A Spartan random process (SRP) is used to estimate the correlation structure of time series and 
to predict (extrapolate) the data values. SRP's are motivated from statistical physics, and they can 
be viewed as Ginzburg-Landau models. The temporal correlations of the SRP are modeled in terms 
of 'interactions' between the field values. Model parameter inference employs the computationally 
fast modified method of moments, which is based on matching sample energy moments with the 
respective stochastic constraints. The parameters thus inferred are then compared with those 
obtained by means of the maximum likelihood method. The performance of the Spartan predictor 
(SP) is investigated using real time series of the quarterly S&P 500 index. SP prediction errors are 
compared with those of the Kolmogorov- Wiener predictor. Two predictors, one of which explicit, 
are derived and used for extrapolation. The performance of the predictors is similarly evaluated. 
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INTRODUCTION 



Time series that carry information about temporal autocorrelations in variables such as 
stock prices, interest rates, etc. have wide applications in finance Such information allows 
predictions of time-series, estimating the associated prediction uncertainty, and performing 
stochastic simulations for reconstructing realizations of the process. For Gaussian time series, 
the temporal structure is estimated from the data by calculating the autocovariance matrix 
or the structure function (variogram). A different approach, typically used in statistical 
physics, focuses on physical interactions embodied in the energy functional. 

Statistical physics plays an increasingly important role in economics, helping to under- 
stand the behavior of complex economical and financial systems (|2j). However, applications 
in the prediction and simulation of such systems are explored less. Recently, a class of Gaus- 
sian random fields, named Spartan Spatial Random Fields (SSRF) [3] was proposed as a 
general framework for spatial modelling. Its main advantages lie in parametric frugality, 
potential for including physical constraints in the probability density function (pdf), and 
efficient model estimation. Herein we define Spartan Random Processes (SRP) in time, 
formulate a Spartan predictor (SP), and investigate its potential for time series prediction. 



SPARTAN RANDOM PROCESSES 



Let us consider a noise-free detrended time-series X\(t) that represents the fluctua- 
tions of an observable with temporal resolution A. For statistical inference from a single 
realization (state), second-order stationarity 4] and ergodicity 0, [s] are often assumed. 
The latter implies that the 'characteristic' scale of the fluctuations be considerably smaller 
than the domain size. In statistical physics, the pdf of stationary Gaussian time series 
can be expressed in terms of an energy functional H[X\(t)}, according to the Gibbs pdf 
/ x [Xa] = Z^ 1 exp {— if [Xa]} , where Z is a normalizing constant (the partition function). 



The SRPs can be defined following the formalism introduced in [3j, restricted to one 
dimensional domains. The fluctuation-gradient-curvature (FGC) SSRF model introduced 

n 

in |3J embodies Gaussian fluctuations and involves three terms that measure the square of 
the magnitude, the gradient and the curvature of the fluctuations. On a ID chain, the FGC 
form of #[A/\] can generally be written as 
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H igc [X x ;0} = ^j dt {iMt)] 2 + Vie[Xx(t)] 2 +e[X x (t)] 2 }, (1) 

where X\(t) and X\(t) denote respectively the first and second time derivatives, and 6 = 
(Voj Vii £■> ^c) is a vector of model parameters: rj is the scale coefficient, 771 the autocovariance 
shape coefficient, £ is the characteristic length, and k c oc A -1 the cutoff frequency. For 
discrete time series sampled at tj = ia, i — 1, . . . , N, a > 0, the energy functional can be 
expressed in terms of local energies S m (t n ), m = 0, 1, 2, as follows: 

1 N 

H igc [X x ; 0] = — V {S (t n ) + m + eS 2 (t n )} , (2) 

where So(t n ) = X\{t n ), Si(t n ) = [(X\(t n + a) — X\(t n ))/a] 2 is the square of the forward- 
difference gradient approximation, and S^n) — [X\{t n + a) + X\{t n — a) — 2X\(t n )} 2 /a 4 
is the square of the discrete approximation of the Laplacian. The ID spectral density is 
given by the following expression 

where K\(k) is the Fourier transform (FT) of the smoothing kernel that imposes the reso- 
lution A. The autocovariance function is obtained from the inverse FT, i.e,. by the integral 

/°° dk ~ 
—G x (k;0)eMjkt). (4) 
-oo 2 7T 

PARAMETER INFERENCE: THE MODIFIED METHOD OF MOMENTS 

The modified method of moments (MMoM) is based on fitting sample constraints 
with corresponding stochastic ones. The former are based on the short-range moments 
S m (t n ), m — 0, 1, 2, appearing in Eq. (T5]). They are evaluated by means of sample averages: 



So(t n ), Si(t n ), and 5 2 (t n ). The respective stochastic constraints can be expressed as follows: 



kc 



^ = ( y„ ) = ^y o (5 ) 

E[ Sl ] = * [(5,(0) - G»] = 2 S / V ll -,™f a L> , (6) 



a 2L v/ v /J ira 2 J 1 + Vl (k£) 2 + (k^ 4 
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£[S 2 ] = - [3G X (0) + G x (2«) - 4G x (a)] = ^J q dk x + ^ ( ^ )2 + ( ^ )4 . (7) 

The stochastic constrains are related to the SRP model and depend on a, r] , rji, £, and 
fc c . For the uniform sampling step a used herein an infinite k c (unlimited-band), will be 
considered, allowing the stochastic constraints Eqs. (EH) to be be expressed in closed form. 
The model parameters are estimated by minimizing the following distance metric {^J 



S~ E[Si 

FCG MODEL BASED SPARTAN PREDICTOR 



S[ E[S 2 



(8) 



Let us assume that T m = {ti,...,tx} is a set of sampling times and X(T m ) = 
{X(ti), . . . , X(tjc)} the set of measurements. The spacing of T m is either uniform (full 
series) or non-uniform (training subset). We assume that T p = {zi, . . . ,Zk} is a set of esti- 
mation points, disjoint from T m , and X(T P ) are the estimated values (temporal predictions). 
T = T m U T p is the full measurement - prediction set with cardinality N = K + k. 

Interpolation: The Spartan predictor (SP) is based on maximizing the conditional 
probability density f x [X(T p )\X(T m )]. Considering the relation f x [X (T p )\X (T m )] = 
fx[X(T)]/fx[X(T m )}, the problem reduces to finding the maximum of f x [X(T)]. We ac- 
complish this by replacing of Hf gc with the estimator Hf gc [X(T);0] = ^X'(T)J x (0)X(T), 
and solving the linear system 



dH {gc [X(T);6]/dX(zi] 



= 0, l = l,...,k. (9) 



Hfgp[X(T);6] involves the sampling points as well as the prediction points, and J x (0) only 
depends on the model parameters, not the data. Neglecting interactions between the pre- 
diction points, the linear predictor can be expressed explicitly by 

X(T P ) = -J2 LeV J,(T L ,T p )X(T L )/J x (T L ,T L ), (10) 

where V is the interaction neighborhood. Herein, the latter extends up to the second- 
nearest neighbor. The numerical complexity of SP involved in solving simultaneously k 
coupled Eqs. <Q is 0(k 3 ). The predictor given by Eq. ffTUj) is explicit. 



Extrapolation: Given the set of measurements X(T m ) = {X(ti), . . . , X(tjv)} a t times 
T m = {ti, . . . , ijv}, we aim to estimate k future values, X(T P ) p = N + 1, . . . , N + k. 
Two approaches are possible: (i) Multipoint (MP) extrapolation, i.e. solving simultane- 
ously the system of equations where I = N + 1, . . . , N + k, or (ii) iterative feed- 
forward (IFF) point-like prediction. The latter is based on the short-range memory property 
fx[X(ti)\X({ti-3,ti-. 2 ,ti-i})] = fx[X (ti)\X ({ti- 2 ,ti-i})}, which allows the following explicit 
predictor, where X(ti) = X(ti) for i — N — 1, N : 

X(t t ) = - [j,_ 2) ,X(t,_ 2 ) + J^XiU^)] I J M I = N + 1, . . . , N + k. (11) 

DESCRIPTION OF DATA AND METHODS OF ANALYSIS 

The time series used in this study consists of 388 quarterly (a = 1/4 year) S&P 500 
index data, recorded in 1900-1996 [7]. The Spartan parameters are estimated by means of 
the MMoM and the maximum likelihood estimation (MLE) method js|, using a training 
sets of 132 points. These are randomly selected from the 388 points to obtain 100 different 



configurations. The optimization uses the Nelder-Mead simplex search algorithm [9J and 
is terminated when both the model parameters and the cost function change between con- 
secutive steps less than e = 10~ 6 . The cost function is given by the negative log-likelihood 
-NLL- function in the case of MLE, and the distance metric -DM-, given by Eq. (jHJ), in the 
MMoM case. Initial guesses for the Spartan model parameters are ^ = a and 77^ = 1. 

To evaluate the prediction performance, for each realization of the training set, the re- 
maining 256 points (validation set) are predicted by both SP (MP and IFF) and Kolmogorov- 



Wiener predictor (KWP) 10]. In both cases, the Spartan covariance model is used. In 
KWP, the search neighborhood includes the entire series. The following statistics are eval- 
uated, where Xi is the real value, Xi is the estimate, and M is the number of validation 
points: (i) mean absolute error (MAE): X]f=i — Xj\/M (ii) mean relative error (MRE): 
J2Zi( X i ~ Xi)/Xi (iii) mean absolute relative error (MARE): J^Zi \( X i ~ X^/MX^ (iv) 



root mean square error (RMSE): yY2i^i \Xi — Xi\ 2 /M, and (v) the linear correlation co- 
efficient (R). The computations are performed in the Matlab® environment on a desktop 
computer with a Pentium 4 CPU at 3 GHz and 1 GB of RAM. 
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FIG. 1: Empirical correlation function (0) 
and optimal Spartan model obtained from 
complete data by means of MLE (★) and 
MMoM (o); for dash-dot line, see text. 

RESULTS 



FIG. 2: Estimates of £ and 771, obtained from 
100 training sets of 132 points by MMoM (o) 
and MLE (★). Note: in the case of MMoM 
some estimates lie outside the axes' range. 



Estimation of Correlations: The empirical correlation function is compared with those 
obtained with the MMoM and the MLE estimators in Fig. [H They all match very well 
near the origin, which is crucial for interpolation. The value of the MMoM distance met- 
ric function is $ = 8.2 x 10 -20 , indicating excellent match of the sample and stochastic 
constraints. The optimization CPU time of MMoM (0.078 s) is 228 times faster than the 
MLE one (17.08 s). Unlike the MLE CPU time that increases nonlinearly with the data, 
the MMoM CPU time is insensitive to the domain size j^]. The difference between the two 
methods in computational time is expected to increase dramatically with the sample size. 

As shown in Fig. [2] the distribution of parameter estimates, obtained from different 
training sets, is scattered and skewed for both MMoM and MLE cases. Some MMoM 
estimates receive extreme values, e.g. r/ > 10 2 , rjx > 10 3 , and £ < 10° (not shown in Fig. 

producing an almost linear decrease of the correlations (e.g., dash-dot line in Fig. [1] 
corresponds to £ = 0.065, rjx = 5 x 10 5 ). Such cases are marked by relatively high $ values, 
indicating poor matching of the sample and stochastic moments, likely due to non-ergodic 
conditions. Indeed, the integral scale [U] is quite large 1(771 = 55.89, £ = 2.72) = 41.4, 
compared to the time series' length. 
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TABLE I: Interpolation errors based on 100 training sets of 132 points. 





MAE 


MARE 


MRE 


RMSE 


R 


MMoM 


0.0717 


0.0763 


-0.0168 


0.1174 


0.9519 


MLE 


0.0731 


0.0771 


-0.0140 


0.1202 


0.9493 



O training set 
x validation set 
• estimated data 




50 100 150 200 250 300 350 

t [1/4 yr] 



FIG. 3: Time series reconstruction by interpolation: training (o), validation (x), estimates (•). 

Interpolation: In the absence of long-range correlations, the short distances are most 
relevant for interpolation. This is evidenced from Table [IJ which compares the interpolation 
errors of the MLE- and MMoM-based estimates: the non-ergodic effect does not introduce 
significant deviations in the MMoM-based estimates. We evaluate the performance of the 
interpolation methods by comparing the observations at the validation points with the pre- 
dictions, using the statistical measures of performance defined above. The model parameters 
estimates are based on the MMoM. In Table [Til the SP prediction errors are compared with 
those of KWP. The comparison distinguishes between areas with different densities of train- 
ing data. Nine categories are defined according to the number of training-point neighbors 
in the vicinity of the validation point. The category implies points from the validation 
set with % nearest and j next-nearest neighbors that belong in the training set. Naturally, 
the categories with more data in their interaction neighborhood display smaller errors than 
those with fewer data. There are no significant differences between the SP and KWP results. 
The linear correlation coefficient values (0.952 for SP and 0.953 for KWP) indicate strong 
correlation between the predictions and the actual data. An example of a reconstructed 
time series from 132 points, obtained by means of SP is shown in Fig. [3j 
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TABLE II: Interpolation errors of SP and KWP predictions based on 100 training configurations. 





(2,2) 


(1,2) 


(0,2) 


(0,1) (0,0) (1,1) 


(1,0) 


(2,0) 


(2,1) 


Total 


Mean Absolute Error 


SP 


0.041 


0.050 


0.063 


0.085 0.126 0.051 


0.053 


0.036 


0.034 


0.0718 


KWP 


0.042 


0.051 


0.063 


0.085 0.123 0.052 


0.053 


0.036 


0.034 


0.0717 


Mean Absolute Relative Error 


SP 


0.038 


0.052 


0.067 


0.097 0.129 0.054 


0.059 


0.039 


0.033 


0.0756 


KWP 


0.040 


0.053 


0.067 


0.097 0.126 0.055 


0.059 


0.039 


0.033 


0.0761 


Mean Relative Error 


SP 


-0.001 


-0.005 


-0.010 


-0.034 -0.032 -0.010 


-0.016 


-0.003 


-0.001 


-0.0168 


KWP 


-0.001 


-0.004 


-0.009 


-0.034 -0.032 -0.010 


-0.015 


-0.002 


0.000 


-0.0160 


Root Mean Square Error 


SP 


0.051 


0.073 


0.089 


0.130 0.174 0.082 


0.084 


0.052 


0.049 


0.1174 


KWP 


0.053 


0.075 


0.089 


0.130 0.165 0.082 


0.084 


0.053 


0.050 


0.1162 



Extrapolation: In Table IIHI we compare SP errors, obtained by the SP-MP and SP-IFF 
methods with those obtained using a second-order autoregressive (AR) model. We randomly 
select from the iV = 388 data a point i (2 < i < N — k) and predict the following k = 3 
values. The errors obtained by the two SP methods are almost identical and comparable 
(slightly smaller) with those obtained from the AR model. At 1-st, 2-nd and 3-rd lag, the 
predicted future values have average relative errors on the order of 5%, 10%, and 15%, 
respectively. 

SUMMARY 

We present a framework for the analysis of time series based on 'pseudo-energy' func- 
tionals that capture the temporal heterogeneity of the observed process. Estimates of the 
process at unmeasured points (predictions) are based on the mode of the joint pdf. The 
Spartan prediction method is tested on financial data (S&P 500 index time series), by both 
interpolation and extrapolation. The SP yields results well comparable to those obtained 
by the standard Kolmogorov- Wiener predictor at generally lower computational cost. 
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TABLE III: Extrapolation error statistics of SP (MP and IFF) and AR model based on 100 
randomly selected triplets i + l,i + 2,i + 3 (2 < i < N — 3). 



Method 




SP-MP 






SP-IFF 






AR 




Lag r 


1 


2 


3 


1 


2 


3 


1 


2 


3 


MAE 


0.0549 


0.0937 


0.1390 


0.0548 


0.0937 


0.1390 


0.0583 


0.1068 


0.1549 


MARE 


0.0495 


0.0833 


0.1383 


0.0495 


0.0832 


0.1384 


0.0533 


0.0996 


0.1633 


MRE 


0.0047 


0.0074 


-0.0156 


0.0048 


0.0075 


-0.0155 


-0.0096 


-0.0211 


-0.0583 


RMSE 


0.0696 


0.1164 


0.1857 


0.0695 


0.1164 


0.1858 


0.0827 


0.1396 


0.2110 



We also present an efficient parameter inference technique - modified method of moments 
(MMoM), which is based on fitting of sample and corresponding stochastic short-range con- 
straints. The main advantage of MMoM is low computational complexity (high speed) and 
independence on the domain size, which makes it suitable for large data sets, difficult to man- 
age by other techniques. Alternatively, the computational ease gives the method potential 
include time-dependent parameters that are continuously estimated, in a "moving window" 
approach, to account for potential non-stationarity. Since the Gaussian assumption is of- 
ten unjustified, we currently focus on formulating a Spartan model capable of representing 
directly non-Gaussian data, without the need for a normalizing transformation. 

This work is partially supported by the Marie Curie Transfer of Knowledge Program, 
Project SPATSTAT 014135. 
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